% Sample Input for Coarse Mesh Diffusion Solver COMEDS2D
% !! This input will correspond to a version which does
%    not use interpolation, generating new responses for
%    each outer iteration !!
% j.roberts 2/23/2010 0.962860682991248, 0.64185, 0.6419
%clear; 
%profile on
name = 'sample_input_v1';

%--------------------------------------------------------------------------
% BASIC PROBLEM INPUT VARIABLES

% number of element types
numtypes = 2;
% number of groups
numg  = 2; 
% number of materials
numm = 2;
% legendre order
order = 2;

%--------------------------------------------------------------------------
%  READ IN DATA -- HERE, JUST DEFINED 
%        D      sigR    nu*sig   xi     scattering
dat = [ 1.500   0.030   0.000    1    0.000  0.000  
        0.420   0.080   0.135    0    0.020  0.000
        2.000   0.040   0.000    0    0.000  0.000
        0.300   0.010   0.000    0    0.040  0.000 ]; 

% dat = [ 1.000   0.50   0.600  1.00  0.000  
%         1.000   0.50   0.800  1.00  0.000 ]; 
abdat=0;   
%--------------------------------------------------------------------------     
%  PLACE THE ELEMENTS
%elements = [ 2 2 2 2
%             2 1 1 2
%             2 1 1 2
%             2 2 2 2];  
elements = [ 1 1 
             1 1 ];  
% 0.90487828

if (max(max(elements))>numtypes)
    disp('numtypes < largest value of elements')
    return
end

bc(1) = 1;
bc(2) = 0;
bc(3) = 1;
bc(4) = 0;

xcm  = [   0   4 ]*5;
xfm  = [     1    ]*10;
ycm  = [   0    4 ]*5;
yfm  = [     1    ]*10;

% compute number in each region
numycm = length(yfm); % number of y regions
numxcm = length(xfm); % number of x regions
mt = zeros(numxcm,numycm,numtypes);

% ELEMENT 1 - material specification
mt(:,:,1)   = [  1 ]; 
% % ELEMENT 2 - material specification
mt(:,:,2)   = [  2 ]; 
%                                       
tol     = [1.e-10 1.e-10]; % tolerance for norms             
maxit   = 10;

swpi = 1;
swnt = 1;
swgm = 0;

input = makeinput( name, numtypes, numg, numm, order, dat, abdat, ...
                   elements, bc, xcm, xfm, ycm, yfm, mt, tol, maxit, ... 
                   swpi, swnt, swgm, 1 ) ; 
              
% return  
%-------------------------------------------------------------------------
% RUN SOLVER and OUTPUT EITHER TO FILE OR SCREEN      

tic
[x1,normk1,it,nrm1] = solverPI_mc(input);
fprintf(1,'  final keff = %12.8f\n',x1(end-1));
toc


tic
input.maxit=1;
[xseed,tmp,tmp,tmp]=solverPI_mc(input);
[x2,normk2,it2,nrm2] = solverNT_mc(input,xseed);
fprintf(1,'  final keff = %12.8f\n',x2(end-1));
toc


